#include "hw4.hh"

void parComputeMorton ( point *x, int N, int b ) {
	double min = x[0].x;
	double h = (x[N-1].x - min) / (double) b;
	int i;

	h += h * 0.0001;

	#pragma omp parallel for
	for (i = 0; i < N; i++)	x[i].m = floor((x[i].x - min) / h); 
}

int parComputeDepth ( double *x, int N ) {
	int i;
	int b;
	double h = 99;
	double d;

	#pragma omp parallel for private(d) 
	for (i = 0; i < N-1; i++) {
		d = x[i+1] - x[i];
		#pragma omp critical
		if ( d < h ) { 
			h = d;
		}
	}

	h *= 0.9999;
	b = pow2roundup((int)((x[N-1] - x[0]) / h));

	return (int)ceil(log2(b));
}

void G ( double *y, double *d, int N, double *x, double *u, int M ) {
	return;
}

TreePtr *setup ( double *y, double *d, int N, 
		void (*G)(double *y, double *d, int N, double *x, double *u, int M)) {
	int i;						// Iterator.
	int depth;					// Depth of tree.
	int buckets;				// Buckets for points.
	int ppt;					// Points per thread.
	int thr;					// Number of running thread
	TreePtr *T;					// Main btree of points.
	TreePtr *TA[numThreads];	// Array of one tree per thread.
	point *y_pt;				// Array of points.
	point null_pt;
	null_pt.m = -1;	null_pt.x = -1;

	parqsort<double>(y, 0, N-1, &cmpDouble, &swapDouble);

//	for (i = 0; i < N; i++) cout << y[i] << " "; cout << endl;
//	depth = (int)ceil(log2(N));
	depth = parComputeDepth( y, N );

	buckets = 1 << depth;
	ppt = ceil( 1.0 * N / numThreads );
	y_pt = new point[N];

	for (i = 0; i < N; i++) {
		y_pt[i].x = y[i];
		y_pt[i].origRank = i;
		y_pt[i].mID = 1;
	}
	
	  // Sort points and compute Morton IDs in parallel
	parComputeMorton(y_pt, N, buckets);

	for (i = 0; i < N; i++) cout << y_pt[i].m << " "; cout << endl;

	  // Initialize trees
	for (i = 0; i < numThreads; i++) {
		TA[i] = new TreePtr( null_pt, depth, buckets );
		TA[i]->Root = TA[i];
	}

	cout << "ppt = " << ppt << ", thr = " << numThreads << endl;
	#pragma omp parallel for schedule(static, ppt) private(thr) 
	for (i = 0; i < N; i++) {
		thr = omp_get_thread_num();
		TA[thr]->Insert( y_pt[i] );
	}

	for (i = 0; i < numThreads; i++) {
		TA[i]->Preorder();	
		cout << endl;
	}
}

int main ( int argc, char **argv ) {
	int N = 20;					// # of points. Default to 20.
	int i;						// Iterator.
	double *y;					// Array of points.
	double *d;

	  // Read in N if given.
	if (argc > 1) N = atoi(argv[1]);
	if (argc > 2) numThreads = atoi(argv[2]);

	  // Allocate space for array of points.
	y = new double[N];

	  // Seed random number generator
	srand48( time(NULL) );

	  // Seed first half of array w/normal dist: u=0.3, v=0.05
	for (i = 0; i < N/2; i++) y[i] = normal(0.3, 0.05);	

	  // Seed second half of array w/normal dist: u=0.6, v=0.01
	for (i = N/2; i < N; i++) y[i] = normal(0.6, 0.01);

	  // Set number of threads
	omp_set_dynamic( 1 );
	omp_set_num_threads( numThreads );

	  // Setup tree
	setup( y, d, N, &G );

	return 0;
}

